#!/usr/bin/env jython

CURRENT_G_TEST, CURRENT_L_TEST, SLOPE_EQUALS = (1, -1, 0)

def tan_test(pre_current, current_vertex, test_vertex):
	"""Returns -1, if the slope of line includes current_vertex and pre_current is greater than the slope of line includes test_vertex and pre_current.
	 Returns 0 if euqal.
	 Returns 1, if the slope of line includes current_vertex and pre_current is less than the slope of line includes test_vertex and pre_current."""
	return cmp((current_vertex[0] - pre_current[0])*(test_vertex[1] - pre_current[1]) - (test_vertex[0] - pre_current[0])*(current_vertex[1] - pre_current[1]), 0)

def _distance(p1, p2):
	"""Returns the squared distance between p1 and p2."""
	dx, dy = p2[0] - p1[0], p2[1] - p1[1]
	return dx * dx + dy * dy

def _keep_left(hull, test_vertex):
	while len(hull) > 1 and tan_test(hull[-2], hull[-1], test_vertex) != CURRENT_G_TEST:
	# The slope of pq must be greater than the slope of rp
		hull.pop()
	if not len(hull) or hull[-1] != test_vertex:
		hull.append(test_vertex)
	return hull

def _graham_scan(points):
	"""Returns points on convex hull of a subhull's points in clockwise order."""
	points = sorted(points)
	l = reduce(_keep_left, points, [])
	u = reduce(_keep_left, reversed(points), [])
	return l.extend(u[i] for i in xrange(1, len(u) - 1)) or l

def _rtangent(hull, pre_current):
	"""Return the index of the point in hull that the right tangent line from pre_current
	to hull touches.
	"""
	l, test_vertex = 0, len(hull)
	l_prev = tan_test(pre_current, hull[0], hull[-1])
	l_next = tan_test(pre_current, hull[0], hull[(l + 1) % test_vertex])
	while l < test_vertex:
		c = (l + test_vertex) / 2
		c_prev = tan_test(pre_current, hull[c], hull[(c - 1) % len(hull)])
		c_next = tan_test(pre_current, hull[c], hull[(c + 1) % len(hull)])
		c_side = tan_test(pre_current, hull[l], hull[c])
		if c_prev != CURRENT_L_TEST and c_next != CURRENT_L_TEST:
			return c
		elif c_side == CURRENT_G_TEST and (l_next == CURRENT_L_TEST or
									  l_prev == l_next) or \
				c_side == CURRENT_L_TEST and c_prev == CURRENT_L_TEST:
			test_vertex = c			   # Tangent touches left chain
		else:
			l = c + 1		   # Tangent touches right chain
			l_prev = -c_next	# Exchange the sides
			l_next = tan_test(pre_current, hull[l], hull[(l + 1) % len(hull)])
	return l

def _min_hull_pt_pair(hulls):
	"""Returns the hull, point index pair that is minimal."""
	h, pre_current = 0, 0
	for i in xrange(len(hulls)):
		# Returns the minmum value of each hulls[i]
		j = min(xrange(len(hulls[i])), key=lambda j: hulls[i][j])
		if hulls[i][j] < hulls[h][pre_current]:
			h, pre_current = i, j
	return (h, pre_current)

def _next_hull_pt_pair(hulls, pair):
	"""
	Returns the (hull, point) index pair of the next point in the convex
	hull.
	"""
	pre_current = hulls[pair[0]][pair[1]]
	next_pt = (pair[0], (pair[1] + 1) % len(hulls[pair[0]]))
	for h in (i for i in xrange(len(hulls)) if i != pair[0]):
		s = _rtangent(hulls[h], pre_current)
		current_vertex, test_vertex = hulls[next_pt[0]][next_pt[1]], hulls[h][s]
		t = tan_test(pre_current, current_vertex, test_vertex)
		if t == CURRENT_L_TEST or t == SLOPE_EQUALS and _distance(pre_current, test_vertex) > _distance(pre_current, current_vertex):
			next_pt = (h, s)
	return next_pt

def Chan_ConvexHull(pts):
	"""Returns the list of points in a clockwise order"""
	for m in (1 << (1 << t) for t in xrange(len(pts))):
	# 1<< t means pow(2,t), where t must be a positive integer
		hulls = [_graham_scan(pts[i:i + m]) for i in xrange(0, len(pts), m)]
		hull = [_min_hull_pt_pair(hulls)]
		for throw_away in xrange(m):
			pre_current = _next_hull_pt_pair(hulls, hull[-1])
			#Test whether it completes a hull
			if pre_current == hull[0]:
				return [hulls[h][i] for h, i in hull],hulls,m
			# If it does not complete a hull, then add this point as a member of the vertices
			hull.append(pre_current)

if __name__ == "__main__":
	import sys
	if '/usr/share/pyshared/' not in sys.path:
		sys.path.append('/usr/share/pyshared/')
	from reportlab.graphics.shapes import Line,Circle,Drawing
	from reportlab.pdfgen.canvas import Canvas
	from reportlab.graphics import renderPDF
	from reportlab.lib import colors
	from reportlab.lib.units import inch
	from random import gauss,random
	from datetime import datetime
	from copy import deepcopy
	if len(sys.argv) < 2:
		numb = 256
	else:
		numb = sys.argv[-1]

	drawing = Drawing(400, 400)
	pts=[(gauss(200,25),gauss(200,25)) for i in xrange(int(numb))]
	output = Canvas("Convex_hull.pdf")
	output.setPageSize((500,500))

	begin = datetime.now()
	result = Chan_ConvexHull(pts)
	chan_s_hull = result[0]
	sub_hulls = result[1]
	H = result[2]
	end = datetime.now()
	delta = end -begin
	#print delta.microseconds,len(chan_s_hull)
	
	if int(numb) <= 256:
		for i in pts:
			drawing.add(Circle(i[0],i[1],0.75,strokeWidth=0))
		
#		for i in sub_hulls:
#			tmp = i+[i[0]]
#			current_color = colors.Color(random(),random(),random(),1)
#			for j in xrange(len(tmp)-1):
#				drawing.add(Line(tmp[j][0],tmp[j][1],tmp[j+1][0],tmp[j+1][1],strokeColor = current_color,strokeWidth = 0.125))
#			renderPDF.draw(drawing,output,50,50)
#			output.showPage()
		
		output.setFont("Helvetica-Bold", 25)
		output.drawCentredString(200, 300, "Algorithm step by step")
		output.showPage()

		renderPDF.draw(drawing,output,50,50)		
		output.showPage()

		output.setFont("Helvetica-Bold", 25)
		output.drawCentredString(300, 300, "Find each subhull")
		output.showPage()

		tmp_draw = Drawing(400,400)
		for i in xrange(0,len(pts),H):
			for j in pts[i:i+H]:
				current_color = colors.Color(random(),random(),random(),1)
				tmp_draw.add(Circle(cx=j[0],cy=j[1],r=0.75,fillColor = current_color, strokeColor=current_color,strokeWidth=0))
		renderPDF.draw(tmp_draw,output,50,50)
		output.showPage()




		for i in xrange(0, len(pts), H):
			tmp_draw = Drawing(400,400)
			for j in pts:
				if j in pts[i:i+H]:
					tmp_draw.add(Circle(cx=j[0],cy=j[1],r=0.75,fillColor = colors.red, strokeColor=colors.red,strokeWidth=0))
				else:
					tmp_draw.add(Circle(cx=j[0],cy=j[1],r=0.25,fillColor = colors.black, strokeColor=colors.black,strokeWidth=0))
			sub_hull = _graham_scan(pts[i:i+H])
			tmp = sub_hull+[sub_hull[0]]
			current_color = colors.Color(random(),random(),random(),1)
			for k in xrange(len(tmp)-1):
				tmp_draw.add(Line(tmp[k][0],tmp[k][1],tmp[k+1][0],tmp[k+1][1],strokeColor = current_color,strokeWidth = 0.125))
				drawing.add(Line(tmp[k][0],tmp[k][1],tmp[k+1][0],tmp[k+1][1],strokeColor = current_color,strokeWidth = 0.125))
			renderPDF.draw(tmp_draw,output,50,50)
			if i== xrange(0,len(pts),H)[-1]:
				output.setFont("Helvetica-Bold", 25)
				output.drawCentredString(200, 400, "All subhulls are found")
				output.showPage()
			else:
				output.showPage()
			if i >0:
				renderPDF.draw(drawing,output,50,50)
				output.showPage()


		
		output.setFont("Helvetica-Bold", 25)
		output.drawCentredString(300, 300, "Find the entire Hull")
		output.showPage()
		
		tmp_draw = deepcopy(drawing)
		for i in pts:
			tmp_draw.add(Circle(i[0],i[1],0.25,strokeWidth=0))
		for i in xrange(0,len(pts),H):
			for j in _graham_scan(pts[i:i+H]):
				tmp_draw.add(Circle(cx=j[0],cy=j[1],r=0.75,fillColor = colors.red, strokeColor=colors.red,strokeWidth=0))
		renderPDF.draw(tmp_draw,output,50,50)
		output.showPage()

		tmp_draw = Drawing(400,400)
		for i in pts:
			tmp_draw.add(Circle(i[0],i[1],0.25,strokeWidth=0))
		for i in xrange(0,len(pts),H):
			for j in _graham_scan(pts[i:i+H]):
				tmp_draw.add(Circle(cx=j[0],cy=j[1],r=0.75,fillColor = colors.red, strokeColor=colors.red,strokeWidth=0))
		renderPDF.draw(tmp_draw,output,50,50)
		output.showPage()


		tmp = chan_s_hull+[chan_s_hull[0]]
		for i in xrange(len(tmp)-1):
			drawing.add(Line(tmp[i][0],tmp[i][1],tmp[i+1][0],tmp[i+1][1],strokeColor=colors.red, strokeWidth=0.75))
			renderPDF.draw(drawing,output,50,50)
			output.showPage()


		output.save()

	def trier(points):
		beg_time = datetime.now()
		c_hull = Chan_ConvexHull(points)
		end = datetime.now()
		delta = end -begin
		return delta.microseconds,len(c_hull[0])
	
#	print "n\tvetices\tms"
#	for number in xrange(1001,21000,2000):
#		spends,h_v = trier([(gauss(200,25),gauss(200,25)) for i in xrange(number)])
#		print i,"\t",h_v,'\t',spends


		#renderPDF.drawToFile(drawing, 'Convex_hull.pdf', "Chan's Convex Hull")
	numb =1000000
	pts=[(gauss(200,25),gauss(200,25)) for i in xrange(numb)]

	begin = datetime.now()
	result = Chan_ConvexHull(pts)
	chan_s_hull = result[0]
	sub_hulls = result[1]
	H = result[2]
	end = datetime.now()
	delta = end -begin
	print delta.microseconds,len(chan_s_hull),numb
	
	begin = datetime.now()
	chan_s_hull = _graham_scan(pts)
	end = datetime.now()
	delta = end -begin
	print delta.microseconds,numb
